---
title: "R4.LandReform"
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/ideo_list.Rdata")
```

# Set covariates
```{r}
covs_list <- c("tx_Hi", "as.factor(DistN)", "ln_PopDens_70", "Cities_10km", "ln_Confs_All", "PrRecH_Adm3", "Isolados1950", "PrWheatAr_1972", "LandIneqRat", "TerrainSlopeIndex", "MedAltitude", "Communist_1975", "MinDistPCP_Map", "log_gnr_cumul_23_Z")
```

```{r}
elections <- c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87")
```

# Land Reform Models
```{r}

baseline_list <- list()  
for (i in 2:7){ # for elections between 1976 and 1987

aggregated_1 <- ideo_list[[1]] %>% select(-election)  # Get vote shares for baseline year (1975)
colnames(aggregated_1)[3] <- paste(colnames(aggregated_1)[3], "1", sep = "_")
 
aggregated_2 <- ideo_list[[i]] %>% select(-election)  # Get vote shares for year of interest
colnames(aggregated_2)[3] <- paste(colnames(aggregated_2)[3], "2", sep = "_")

# First Difference
aggregated_fd <- aggregated_2 %>% inner_join(aggregated_1) %>% rowwise() %>% # Subtract 1975 vote shares from vote shares of year of interest
  mutate(FD = votes_2 - votes_1) %>%
  select(DICOFRE_Z, ideology2, FD) %>%
  pivot_wider(c(DICOFRE_Z), names_from = "ideology2", values_from = "FD")

data_sample_g <- data_final %>% left_join(aggregated_fd) %>%  # Join vote share first difference and cumulative GNR events up to election year to main data set
  left_join(gnr_panel %>% filter(elec_number == i)) 

output <- list() # Subset for (3) groups (land reform intensity)
for (g in 0:2){ 
 
output[[g+1]] <- data.frame(election = elections[i], group = g,  model = colnames(aggregated_fd)[-c(1)], tx_Hi = NA, se = NA, n = NA) %>% filter(model != "NA")

      # Run models for each ideology category (4)
      for (r in 1:nrow(output[[g+1]])){ # r<-1
      m.1 <- lm_robust(as.formula(paste(output[[g+1]]$model[r], paste(c(covs_list), collapse= " + "), sep = " ~ ")),
                data = data_sample_g %>% filter(exp_cat2 == g), weight = weights, cluster = subclass)
      
           output[[g+1]][r, "tx_Hi"] <- m.1$coefficients["tx_Hi"]
           output[[g+1]][r, "se"] <- m.1$std.error["tx_Hi"]
           output[[g+1]][r, "n"] <- length(m.1$fitted.values)
      }

}

output_list <- do.call(rbind, output)
baseline_list[[i]] <- output_list

}

# Plot Results

baseline_list_all <- do.call(rbind, baseline_list) %>% 
  mutate(model = ifelse(model == "Center", "Center-Left", model), 
         model = ifelse(model == "Left", "Far Left", model))

years <- data.frame(election = c("AR76", "AR79", "AR80", "AR83", "AR85", "AR87"), 
           year = as.factor(c(1976, 1979, 1980, 1983, 1985, 1987)))
baseline_list_all <- baseline_list_all %>% left_join(years)
baseline_list_all$model <- factor(baseline_list_all$model, level = c("Communist", "Far Left", "Center-Left", "Right"))


ggplot(baseline_list_all, aes(x = year, col = as.factor(group))) +
  geom_point(aes(y = tx_Hi, shape = as.factor(group)), position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.975, n)*se), ymax=tx_Hi + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.95, n)*se), ymax=tx_Hi + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  geom_hline(aes(yintercept = 0), col = "red") +
  facet_wrap(~model, scales = "free_y") +
  xlab("Election Year") + ylab("Treatment Effect") +
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(16, 17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))

ggsave(plot = last_plot(),
       "Figures/fg7_MainReform.pdf",
       width = 9, height = 5)
  
```
